Spin gap and superconductivity in the three-dimensional attractive Hubbard model 
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We study the phase diagram for the attractive (i.e., negative-?/) Hubbard model on a simple cubic 
lattice, through Monte Carlo simulations. We obtain the critical temperature, T c , for superconduc- 
tivity from a finite-size scaling analysis of the data for the pairing correlations. For fixed on-site 
attraction, U, T c displays a maximum near the filling factor 0.9, roughly independent of U. For 
fixed filling we estimate the crossover temperature T X (Z7), separating the normal states: metallic 
and spin-gap. There is also a critical value U p for pair formation, the magnitude of which seems to 
be independent of doping. The relevance of these results to the high-T c oxides is discussed. 
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One of the most striking normal state properties of 
the high-T c cuprate superconductors is the behavior of 
the uniform magnetic spin susceptibility, Xs, as the 
temperature is lowered: Instead of being temperature- 
independent as in conventional Fermi liquids, Xs starts 
to decrease well above the critical temperature T c , as ev- 
idenced by NMR Knight shifts and relaxation rates 1 and 
by direct susceptibility measurements 2 . This suppression 
of Xs has been associated with the opening of a spin gap 
at a crossover temperature, T x , above T c 3 ~ 5 . As the tem- 
perature is decreased further, and within a certain range 
of doping, the material becomes superconductor. In the 
search for a mechanism responsible for superconductivity 
in these materials, it is therefore instructive to study sim- 
plified models displaying the essential features observed, 
such as the precursor spin-gap phase in the normal state. 
The attractive Hubbard model (i.e., negative on-site cou- 
pling U) is believed to display these features 6-8 ' 4 . Early 
mean-field calculations 6 indicate that local singlet pairs 
are formed at high temperatures, and that these inco- 
herent pairs condense into a charged superfluid at T c . In 
two dimensions, Randeria et al. 4 have provided numerical 
evidence to show that the uniform susceptibility of the 
model is suppressed above the superconducting temper- 
ature and that it is proportional to the NMR relaxation 
rate. Due to their layered structure, one should expect 
some properties of the cuprates to interpolate between 
the two- and three-dimensional models. Local fermion 
pairs may be formed in narrow-band systems due to a lo- 
cal attractive short-ranged effective interaction and have 
also been invoked to explain a variety of other phenom- 
ena; see Ref. 6 for a list of references. More recently, 
a model for the Cu02 planes with interacting carrier 
and insulating bands and repulsive interactions has been 
mapped onto the attractive Hubbard model 9 . In view 
of all this, a systematic study of the attractive Hubbard 
model in three dimensions is in order. In particular, there 
are many aspects such as the behavior of T c and T x with 
both U and the occupation away from half-filling, that 
are known at most qualitatively. With this in mind, here 
we address these questions through Monte Carlo simula- 



tions. 

The Hubbard Hamiltonian can be written as 

n = - t Y, ( C «W + H - c -) + u E( n «t - \) ( n n ~ \) 

-/i^n iCT , (1) 

where the sums run over sites of a simple-cubic lattice, 
denotes nearest neighbor sites, H.c. stands for Her- 
mitian conjugate, and c irT (ci a ) creates (annihilates) a 
fermion at site i with spin c; U < is the attractive 
on-site interaction and /i is the chemical potential con- 
trolling the band-filling. Since the simple-cubic lattice is 
bipartite, the band is half-filled when the Hamiltonian (1) 
displays particle-hole symmetry, or /i = 0. In this case, 
superconducting correlations in the attractive model are 
equivalent to planar magnetic correlations in the repul- 
sive model 6 . The strong-coupling limit of (1) can be ob- 
tained through perturbation theory in the space of dou- 
bly occupied states and is equivalent 10 ' 11 to a Heisenberg 
model in a transverse field proportional to /i. 

Here we use a grand-canonical Quantum Monte Carlo 
simulation; see Refs. 12-15 for details. The imagi- 
nary time is discretized through the introduction of M 
"time" slices separated by an interval At such that 
P = At M. One should stress that the simulation for 
the attractive Hubbard model is free from "minus sign" 
problems 11 ' 14 ' 15 . We calculate quantities such as the 
equal-time q = local (or s-wave) pairing correlation 
function, 

P S (T, L) = (A' 1 ' A + AA' 1 ') , (2) 

where T is the temperature, L is the linear lattice size, 
and 

and the uniform magnetic susceptibility 
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The dependence of the pairing correlation function 
with the system size can be extracted through finite- 
size scaling (FSS) arguments 16 . For an infinite three- 
dimensional system one expects a superconducting tran- 
sition within the XY -model universality class, with pair- 
ing correlations decaying algebrically at the critical tem- 
perature T c . For a finite system of size L, one can as- 
sume the following FSS ansatz for its associated uniform 
Fourier transform 11 : 



P S (T,L) = L 2 ~ v F (L I , 



(6) 



where £ is the correlation length for the infinite system, 
and F(z) is a scaling function such that F(z) — ^ const 
when L <C £; in three dimensions 17 , r\ ~ 0. At T c , £ = 
oo, so that L~ 2 P a (T c , L) is a constant independent of 
lattice size. By plotting L~ 2 P S (T,L) as a function of T 
for systems of different sizes, an estimate of T c can be 
obtained as the temperature where two successive curves 
intercept 18 . 




FIG. 1. Size-scaled q = Fourier transform of the s-wave 
pairing correlation function as a function of the inverse tem- 
perature, for lattices with L = 3.17 (squares), L = 4 (circles), 
and L = 6 (triangles), with U = —6 and at half-filling. Ex- 
cept where shown, the error bars are smaller than the data 
points, and the solid lines are guides to the eye. 



The clusters used here have N s = L x x L y x L z sites, 
with periodic boundary conditions; that is, each site is 
connected with its six nearest neighbors through a hop- 
ping term. The simulations were performed on Sun and 
IBM RISC-6000/525 workstations; a single datum point 
involves between 500 and 4000 MC sweeps over all time 
slices and we took At = 0.125. In a grand-canonical 
simulation, for each temperature the chemical potential 
is adjusted to obtain the desired occupation, p = (n). 



Since we are interested in several values of both U and 
p, we had to restrict ourselves to small systems due to 
our limited computer capabilities. From now on, energies 
will be expressed in units where the hopping t = 1, and 
we also set the Boltzmann constant k& = 1. 
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FIG. 2. Critical temperature T c as a function of the mag- 
nitude of the on-site coupling constant for half-filled band; 
below T c the system displays both superconductivity (SUC) 
and charge (CDW) ordering. The results from this work, us- 
ing L13 (solid circles) and L2 (solid squares) are compared 
with those obtained from the repulsive model: Monte Carlo 
simulations (open triangles; Ref. 19) and variational calcula- 
tions (solid and dashed lines; Ref. 20). The dotted line is the 
crossover temperature T x (see text). 

We considered lattices with 4x4x2 and 4x4x4 
sites; but in order to assess possible finite-size effects 
we have also performed a few runs on a 6 x 6 x 6 lat- 
tice. For the L x X L y X L z lattices, one may think of 
several definitions 16 for a mean linear size L, such as 



2.83, 
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(L x L y L z 



L 2 = 3/(L" 
~ 3.17. Figure 



U = ^J3/(L X 2 + Ly 2 + L- Z 

Ly 1 + L~ v ) = 3, and L 3 = 
1 shows the size-scaled pairing correlation, Eq. (2), as a 
function of the inverse temperature, for U = —6 at half- 
filling and for three different lattice sizes; the data for 
the 4x4x2 lattice are plotted assuming L = Li = 3. 
One can define the inverse transition temperature, j3 c , as 
the value where data points for two different-sized sys- 
tems (L,L') superimpose 18 . This implies T c ~ 0.25 and 
T c ~ 0.23 for (4,3) and (6,4) scalings, respectively; us- 
ing L = L 3 = 3.17 for the smaller lattice, one would get 
T c ~ 0.3 from the (4,3) scaling. The definition L = 3 for 
the smaller system then provides estimates for T c that are 
closer to the more reliable scaling (6,4) than L\ or L3. 
Idealy a definite trend would only be detectable through 
the use of systems larger than L = 6, which would be- 
come prohibitively expensive in terms of computer time. 
Taking into account the error bars for the data points, the 
above criterion implies a typical error A(3 C ~ ±1. This 
procedure is repeated for other values of U, to obtain 
T r (U) at half-filling. In Fig. 2, the solid symbols repre- 
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sent the critical temperatures obtained from a (4,L) scal- 
ing, both with L = L2 (squares) and with L = L3 (cir- 
cles). One should have in mind that at half- filling the or- 
dered phase below T C (U) corresponds to both supercon- 
ductivity and charge ordering, since the order parame- 
ter displays full three-dimensional rotational symmetry 6 . 
Also, the attractive model at half-filling is equivalent to 
the repulsive model, with the superconducting and charge 
order parameters becoming the planar (XY) and axial (Z) 
staggered magnetizations, respectively. In view of this, 
in Fig. 2 we compare results for T c from the attractive 
model with the Neel temperature Tjv for the repulsive 
model obtained from very extensive simulations 19 (within 
a different extrapolation to estimate TV; open triangles), 
and from a Gut z wilier- type variational calculation 20 (the 
solid line is the "bare" result (Tjv (£/)), and the dashed 
line is an adjustment ((3.83/6.0) x Tjv) to fit the mean- 
field result to that of high-temperature series expan- 
sions for the Heisenberg model, according to which 17 
Tjv ~ 3.83t 2 /|f7|). The estimates for T c using L\ lie be- 
low the "normalized" Tjf(U) which is probably a lower 
bound; from now on all quoted estimates for T c will be 
based on L^. For weak couplings (i.e., \U\ << t), the 
system is in a BCS-like regime; the difference with re- 
spect to the standard BCS theory being due to the fact 
that quasi-particles with any wavevector can pair, not 
only those close to the Fermi level. Accordingly, the crit- 
ical temperature is still exponentially small 6 , but with 
a different energy scale: T c ~ W exp(— W/|C/|) , where 
W = 3t is one half of the bandwidth. 
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FIG. 3. Uniform susceptibility as a function of temperature 
at half filling for a simple-cubic lattice with L = 4. The 
symbols refer to the values of U shown, and error bars are 
smaller than the data points. 

In Fig. 3 we present the uniform susceptibility as a 
function of temperature for the L = 4 lattice at half- 
filling, and for several values of U. The solid line is 
the non-interacting result, xi°\ f° r t ne same lattice size; 
its divergence as T — t is due to the finite-size of the 
system, since Xs approaches the Pauli behavior if the 
L — ► 00 limit is taken before T — ► 0. For U = —0.5 



and U = —1, the behavior of Xs is similar to that of the 
non-interacting (metallic) case. In contrast, for U < —2, 
the uniform susceptibility is suppressed below some tem- 
perature T X (U). This can be understood in the strong- 
coupling regime by noticing that local pairs are being 
formed and that spin excitations necessarily imply pair 
breaking with an energy cost (gap) of order \U\. The 
formation of local pairs, and the associated spin gap, 
should be reflected in the magnetic properties: bound 
singlet pairs must have smaller response to a uniform 
field than isolated fermions. At intermediate couplings, 
this behavior can be explained along similar lines, in 
terms of pairing correlations 4 . Therefore, T X (U) rep- 
resents a crossover temperature separating two normal- 
state regions: metallic and spin-gap. In Ref. 4, this 
crossover temperature was defined as the one at which 
Xs deviates from a renormalized Random Phase Approx- 
imation. Here we choose a different definition, which 
follows closely the experimental criterion based on NMR 
relaxation measurements, namely as the temperature at 
which Xs is maximum; see e.g., the discussion in Ref. 
5. The crossover line obtained this way is displayed in 
Fig. 2. We have compared the data in Fig. 3 with some 
obtained for the L = 6 lattice, and found no significant 
finite-size effects. Nevertheless, in view of the arbitrari- 
ness of this definition, the crossover line obtained is only 
a crude estimate, and should be taken with care. 
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FIG. 4. Critical temperature as a function of the on-site 
coupling —U/t, for p = 0.8 (solid squares); the solid line is a 
guide to the eye. The dotted line is the crossover temperature 
T x (see text). 

As the occupation is varied, the behavior does not 
change drastically. For instance, in Fig. 4 we show 
both T x and T c as functions of — U, for an occupation 
p = 0.8. Away from half-filling the order parameter 
is taio-dimensional, corresponding to superconductivity 
alone; i.e., charge ordering is lost. While T x is roughly 
the same (within the range of U examined here) as for 
p = 1, T c is slightly higher than for the half-filled case. 
A plot of p(fJi) for the non-interacting L = 4 system at 
zero temperature displays several plateaux; in particular 
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there is one at p = 0.6875, corresponding to 44 fermions 
in the system. These plateaux are still present in the 
interacting case, and are rounded at finite temperatures. 
This is a finite-size effect that should disappear in the 
thermodynamic limit, but nonetheless affect the data for 
p = 0.7 in these small systems: the uniform susceptibil- 
ity is strongly suppressed for any U < 0. The data for 
p = 0.6 are free from these effects. 
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FIG. 5. Critical temperature as a function of the band fill- 
ing, for the values of U labeling the curves. The solid lines 
are drawn to guide the eye. 

In Fig. 5 we show the dependence of T c with p for 
several values of U. The data are consistent with T c 
displaying a maximum around p ~ 0.9. For the range of 
U studied, T c increases with \U\ for fixed occupation, but 
should eventually decrease in the strong-coupling limit. 
In the dilute limit (i.e., p —¥ 0), T c should approach 
zero, for any U. But Fig. 5 indicates that the range of 
fillings for which finite-temperature superconductivity is 
effective increases with \U\. 

It is interesting to note 6 that pairs are not formed 
for any U, but below a critical value U p , a precise es- 
timate of which would require further extensive simula- 
tions. Nevertheless, by inspection of Fig. 3 one can say 
that Up £ [— 2,— f], since Xs is suppressed for U = —2, 
but not for U = — 1 ; this value is quite smaller than the 
value Up ~ 7.8 for any p, predicted within a low-density 
approximation 6 as the binding energy of the pair. For 
p=0.9, 0.8 and 0.6, U p lies in the same interval, sug- 
gesting that Up may be insensitive to the occupation. 
The crossover temperature for U = —2 and p = 0.6 is 
about 30% smaller than for the other fillings, while for 
U € [—6, —4] it seems to be less dependent on the occu- 
pation. 

In conclusion, we have obtained the phase diagram in 
the temperature-coupling constant-occupation space for 
the attractive (i.e., negative-?/) Hubbard model on a sim- 
ple cubic lattice. For fixed U, the critical temperature for 
superconductivity displays a maximum at the occupa- 
tion p ~ 0.9. For fixed occupation there are two regimes: 
weak coupling (|J7| <C t), where superconductivity sets in 



at very low temperatures, from a normal metallic state; 
and intermediate- to strong-couplings, where supercon- 
ductivity sets in from a spin-gap phase at higher temper- 
atures. The changeover from a normal metal to a spin- 
gap phase occurs at a crossover temperature, which does 
not seem to be very sensitive to the occupation in the 
range [0.6,1.0], at least for U < —4. regime. We have also 
established that the critical value of \U\ for pair forma- 
tion lies in the interval [—2, —1], for all fillings examined. 
With respect to the cuprates, the existence and origin of 
the superconducting gap has not been fully settled yet. 
If the spin gap, which opens above (and not at) T c in 
underdoped samples, is the only one present, then de- 
scribing superconductivity as arising from the condensa- 
tion of pre-formed pairs as in the model considered here 
is quite appealing. In this respect, we should comment 
on a recent suggestion 21 that the spin gap in the attrac- 
tive model may be irrelevant to the cuprates, as the ob- 
served suppression of the uniform susceptibility would be 
due solely to antiferromagnetic fluctuations. It may be 
possible to reconcile both views if one considers a Hub- 
bard model with on-site repulsion and nearest-neighbor 
attraction. In this case, the superconducting phase is 
near a spin-density wave (SDW) instability 6 , and SDW 
fluctuations could influence the magnetic response as sug- 
gested. Moreover, the superconducting order parame- 
ter in that region has been predicted 6 to have d-wave 
symmetry, in agreement with penetration depth 22 , and 
photoemission 23 studies. Work is in progress to test these 
ideas. 
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